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We present the first results of next-to-leading order QCD corrections to three jet heavy quark 
production at LEP including mass effects. Among other applications, this calculation can be used 
to extract the bottom quark mass from LEP data, and therefore to test the running of masses as 
predicted by QCD. 
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The decay width of the Z gauge boson into three jets has already been computed at the leading order (LO) including 



ihe decay width 01 the Z gauge boson into three jets has already been computed at the leading order (LU j including 
complete quark mass effects Jl|-g[ where it has been shown that mass effects could be as large as 1% to 6%, depending 
on the value of the mass and the jet resolution parameter y c . In fact, these effects had already been seen in the 
experimental tests of the flavor independence of the strong coupling constant [0~0L I n y i ew of that we proposed ||, 
together with the DELPHI collaboration ||, the possibility of using the ratio 
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qq ' as a means to extract the bottom quark mass from LEP data. In this equation T^Jy c ) /T q is the three-jet fraction of 
Z decays into the quark q and y c is the jet resolution parameter. 

Since the measurement of R^ d is done far away from the threshold of b quark production, it will allow, for the first 
time, to test the running of a quark mass as predicted by QCD. However, in [|| we also discussed that the leading 
order calculation does not distinguish among the different definitions of the quark mass, perturbative pole mass, Mb, 
running mass at Mb, or running mass at mz- Therefore in order to correctly take into account mass effects it is 
necessary to perform a complete next -to-leading order (NLO) calculation of three jet ratios including quark masses 

^ MM 

♦pH. In this letter we sketch the main points of this calculation, leaving the details of the complete calculation for other 
publications (l^Jl^], and we present the results that have been used by the DELPHI collaboration to measure the 
D , running mass of the bottom quark at /1 = mz P^JrI . 

• In the last years the most popular definitions of jets are based on the so-called jet clustering algorithms. These 
algorithms can be applied at the parton level in the theoretical calculations and also to the bunch of real particles 
observed at experiment. In the jet-clustering algorithms jets are defined as follows: starting from a bunch of particles 
with momenta pi one computes, for example, a quantity like j/y- = 2 imu(Ef, Ej)/s (1 — cos 9ij) for all pairs (i, j) 
of particles. Then one takes the minimum of all and if it satisfies that it is smaller than a given quantity y c (the 
resolution parameter, y-cut) the two particles which define this are regarded as belonging to the same jet, therefore, 
they are recombined into a new pseudoparticle by defining the four-momentum of the pseudoparticle according to 
some rule, for example, pk = Pi + Pj- After this first step one has a bunch of pseudoparticles and the algorithm can 
be applied again and again until all the pseudoparticles satisfy yy > y c . The number of pseudoparticles found in 
the end is the number of jets in the event. This procedure leads automatically to IR finite quantities because one 
excludes the regions of phase space that cause trouble. It has been shown that, for some of the algorithms, the passage 
from partons to hadrons (hadronization) does not change much the behavior of the observables |fT7f , thus allowing to 
compare theoretical predictions with experimental results. 

Although we have studied the four jet-clustering algorithms discussed in [f^ fn^JT^Jlsfl , here we will present results 
only for the DURHAM algorithm [ p^pCj ] which is the one we just have defined and seems to be the one that behaves 
better for most of the observables. 

The decay width of the Z boson into three jets containing the bottom quark mass can be written as follows 0] 
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= m z / as{mz) (g 2 v H v (yc,r b )+g 2 A H A (y c ,r b )) , (2) 

where g is the SU(2) gauge coupling constant, cw — cos#vy an< 4 s\v = sm8\y are the cosine and the sine of the weak 
mixing angle, gy = — 1 + 4/3s^ and gA = 1 are the vector and axial coupling of the Z boson to the bottom quark, 
and Hv(y c ,rb) and HA{yc,r b ) contain all the dependences in the jet resolution parameter, y c , and the quark mass, 
r b = (M b /mz) 2 , for the vector and axial parts in the different algorithms. These functions can be expanded in a s , 
and, if the leading dependence on the quark mass is factorized out we have, 



H V {A) (ifc, n) = (y c ) + (y c ) + 

r b (B { °} A) (y c ,r b ) + ^B { v %(y c ,r b ))+--- . (3) 

Here, A^°'(y c ) is the tree level contribution in the case of massless quarks and it is known for the different algorithms 
in analytic form. It is exactly the same function for the vector and the axial parts. The functions By} A ^(y c , r b ) take 
into account tree level mass effects once the leading dependence in r b has been factorized out. They were calculated 
numerically in @] for the different algorithms and results were presented in the form of fits to the numerical results. 
^ (j/c) gives the QCD next -to-leading order correction in the case of massless quarks and to good approximation it 
is the same for the vector and the axial parts [Because of the triangle anomaly there are one-loop triangle diagrams 
contributing to Z —* qqg with the top and the bottom quarks running in the loop. Since rut ^ m b the anomaly 
cancellation is not complete. These diagrams contribute to the axial part even for m q = and lead to a deviation 

from Ay\y c ) = A A \y c ) This deviation is, however, small 0] and we are not going to consider its effect here]. 



This function is also known for the different algorithms |L7 18 1. Finally, the functions By^(y ci r b ) were completely 



unknown and contain the next-to-leading order corrections depending on the quark mass. In the next section we 
present our results in the form of fits to some combinations of the relevant functions in the case of the DURHAM 
algorithm, which is the one that gives smaller radiative corrections, and we postpone the presentation of results for 
the different functions and algorithms for another publication Note that the way we write Hv{A)(ya r b) in eq. (||) 

is not an expansion for small r b . We keep the exact dependence on r b in the functions By^ A ^(y c , r b ) and By) A ~. (y c , r b ). 
Factoring out r b makes it easier to analyze the massless limit and the dependence of the results on r b in the region of 
interest. This means that our results can also be adapted, by including the photon exchange, to compute the e + e~ 
cross section into three jets out of the Z peak at lower energies or at higher energies in top quark production. 

At the NLO we have contributions to the three-jet cross section from three and four parton final states. One 
loop three-parton amplitudes are both IR and UV divergent. Therefore, some regularization procedure is needed. 
We use dimensional regularization for both IR and UV divergences because dimensional regularization preserves the 
QCD Ward identities. The UV divergences, however, can be easily removed by renormalization since the appro- 
priate counterterms are very well known. The three-jet cross section is obtained by integrating both contributions, 
renormalized three-parton and four-parton amplitudes, in the thrcc-jct phase space region defined by the different 
jet clustering algorithms. This quantity is infrared finite and well defined, although the three and the four parton 
transition amplitudes independently contain infrared singularities. 

The three-parton transition amplitudes can be expressed in terms of a few scalar one-loop integrals. The result 
contains poles in e — (4 — D)/2, where D is the number of space-time dimensions. Some of the poles come from UV 
divergences and the other come from IR divergences. After UV renormalization we obtain analytical expressions for 
the terms proportional to the infrared poles and for the finite contributions. The infrared poles will cancel against 
the four-parton contributions. The finite contributions are integrated numerically in the three-jet region. 

The four-parton transition amplitudes are split into a soft and collinear part in the three-jet region and a hard 
contribution. The soft terms are integrated analytically in arbitrary D dimensions in the region of the phase space 
containing the infrared singularities. We obtain analytical expressions for the infrared behavior of the four-parton 
transition amplitudes and we show how these infrared terms cancel exactly the infrared singularities of the three- 
parton contributions. The hard terms are calculated in D = 4 dimensions. The remaining phase space integrations, 
giving rise to finite contributions, areperformed numerically. 

Following Ellis, Ross and Terrano [B2j we have classified both, three-parton and four-parton transition probabilities, 
according to their color factors. It is clear that the cancellation of IR divergences between three-parton and four-parton 
processes can only occur inside groups of diagrams with the same color factor. The cancellation of IR divergences 
can be seen more clearly by representing the different amplitudes as the different cuts one can perform in the three- 
loop bubble diagrams contributing to the Z-boson selfenergy. After summing up the three-parton and four-parton 
contributions to the three-jet decay width of the Z boson we obtain the functions Hy(y c , r b ) and H A (y c , i"b) m ec l- (§) 
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at order a s . Details of the calculation, cancellation of divergences and results for the relevant functions will be 
presented elsewhere @Q. Since a large part of the calculation has been done numerically, it is important to have 
some checks of it. We have performed the following tests: 

• We have checked our four parton amplitudes in the massless limit against the amplitudes presented by Ellis- 
Ross- Terrano (ERT) pi| . The three-parton amplitudes for massive quarks cannot be compared directly with 
the massless results because collinear poles in e, in the massless case, appear as logarithms of the quark mass 
when the mass is taken into account. 

• The four parton transition amplitudes have also been checked in the case of massive quarks, in four dimensions, 
by comparing their contribution to four jet processes to the known results Q 

• Of course we have checked that, also in the massive case, all the IR divergences cancel between three parton 
and four parton contributions p3l-E5| . 



• To check the performance of the numerical programs and the overall approach we have applied our method to 
the massless amplitudes of ERT and obtained the known results for the functions AW . 

• We have checked, independently for each of the groups of diagrams with different color factors, that the final 
result obtained with massive quarks reduces to the massless result in the limit of very small masses. 

The last test is the main check of our calculation. We have calculated the functions Hviucfb) and HA{y C if b ) for 
several small values of r b , in the range M b ^1 — 5 GeV, and then we have extrapolated the results for r b — > 0. 
In that limit we recover the values for the function A^\y c ) in the different algorithms considered and the different 
groups of diagrams. This check is not trivial at all since the structure of IR divergences for massive quarks is quite 
different from the case of massless quarks: for massive quarks collinear divergences are regulated by the quark mass, 
and therefore some of the poles in e that appear in the massless case are softened by logn,. Although these checks do 
not constitute a complete test of the massive corrections, we think that all together give some kind of confidence in 
the final result for massive quarks. 

Now, to obtain we can substitute eq. (|^) in eg. (p), this into eq. ([[]) and use the value of T b which is also well 
known (for mass effects at order a s see for instance P,p6 ,p7{ ) . Putting everything together and expanding in a s and 
r b we can write 

Rf=l + r b (b {y c ,r h ) + ^b 1 (y c ,r b j) , (4) 



where the functions bo and b\ are an average of the vector and axial parts, weighted by cy = gy/{gy + g\) and 
ca = 9a/{9v + 9a) respectively, and can be written in terms of the different functions introduced before, eq. (0), 

§0 . . 

It is important to note that because the particular normalization we have used in the definition of R$ , which is 

manifested in the final dependence on cy and ca, most of the electroweak corrections cancel. Those are about 1% 
p8f in total rates while in R^ d are below 0.05%. Therefore for our estimates it is enough to consider tree level values 
of gv and gA- The same argument applies for the passage from decay widths to cross sections. Contributions from 
photon exchange are small at LEP and can be absorbed in a redefinition of g v and g\ f2(| . They will add a small 
correction to our observable. 

Although intermediate calculations have been performed using the pole mass we can also re-express our results in 
terms of the running quark mass by using the known perturbative expression 



Mi = fhiQj) 



tt \ 3 ur 



(5) 



The connection between pole and running masses is known up to order a 2 , however consistency of our pure perturbative 
calculation requires we use only the expression above. We obtain 



R b 3 d = l + f b (n) [b (y c> r b ) + ^- {b 1 (y c ,r b )-2b (y c ,r b )log^-)) . (6) 
Where f b (fi) — m\ (fi) /m 2 z and 

b 1 (y c ,r b ) = b 1 (y c ,r b ) + b Q (y c ,r b )[8/3-2log(r b )} . (7) 



3 



f&(/x) can be expressed in terms of the running mass of the b quark at /i = mz by using the renormalization group. 
At the order we are working 

nW = n{mz) 7-T- with a a (p) = — — -. r-— 8 

V <MMJ / l + a s (m z )p t 

and t = log(^ 2 /m|)/(47r), /3 = 11 - 2N f /3, N f = 5 and 7o = 2. 

At the perturbative level eq. (Q) and eq. (^|) are equivalent. However, since they neglect different higher order terms 
they lead to different answers. Since the experiment is performed at high energies (the relevant scales are mz and 
Tnzy/yc) one would think that the expression in terms of the running mass is more appropriate because the running 
mass is a true short distance parameter while the pole mass contains in it all the complicated physics at scales \i ~ 
Moreover, by using the expression in terms of the running mass we can vary the scale in order to estimate the error 
due to the neglect of higher order corrections. In any case, if one would push the result of eq. (p) up to scales as low 
as fx = 5 GeV one would get something closer to the pole mass result. Therefore, we use eq. (||) and vary the scale 
in an appropriate range to obtain an estimate of the error in the calculation. 

The function bo(y Cl r^) gives the mass corrections at leading order. As shown in || it depends very mildly on the 
quark mass in the region of interest ( Mb ^3 — 5 GeV). Therefore it is appropriate to present our results for bo(y c , r^) 
as a fit in only y c : ba(y c ,rb) — X)n=o ^o™' > l°S™2/c- For the DURHAM algorithm, in the range 0.01 < y c < 0.10 and 
3 GeV <M b <5 GeV and using sin 2 6 W = 0.2315 we obtain k ( ° ] = -10.521 , k { 1] = -4.4352 , = -1.6629 . 

The function bi(y c ,rb) is the main result of this paper. It gives the NLO massive corrections to our observable. It 
is important to note that b\ contains significant logarithmic corrections depending on the quark mass. They come 
from different sources: first, the NLO corrections written in terms of the pole mass now contain some residual mass 
dependence, second the normalization to the total rate induces some additional logarithmic dependences, and finally 
the passage from the pole mass to the running mass adds also an additional logarithmic dependence. Therefore, we 
choose to include explicitly this dependence on the quark mass in our fits to the function bi(y Cl rt,). For our purposes a 

fit of the form b\ (y c , r^) = fcj ' + k[^ log(y c ) + k$ log(r&) is good enough. The coefficients we obtain for the DURHAM 

algorithm and in the ranges we just mentioned for y c and r b are: fcj 0) = 297.92 , k[ 1] = 59.358 , k$ = 46.238 . 

In fig. |l| we present R\ d for // = mz (dashed), fi — 30 GeV (dashed-dotted) and fi = 10 GeV (dotted) for 
mb(mz) = 3 GeV and a s (mz) = 0.118. For comparison we also present the LO results for Mb = 5 GeV (lower solid 
line) which is, roughly, the value of the pole mass obtained at low energies and fhb{mz) = 3 GeV (upper solid line) 
which is, roughly, the value one obtains for the running mass at the mz scale by using the renormalization group 
pf| . Note that choosing a low value for /i makes the result closer to the LO result written in terms of the pole mass, 
while choosing a large /i makes the result approach to the LO result written in terms of the running mass at the 
m z scale. If R\ d is measured to good accuracy one could use eq. (^|) and eq. (||) to extract rhb(mz)- However, the 
extracted result will depend on /i. For illustration, in fig. ^ we represent, as a function of fi, the value one would 
obtain for fhb(mz) if R^exp — 0-96 for y c — 0.02. What is the best scale one should choose in three-jet quantities is a 
long standing discussion. One would think that if the energy is equally distributed among the three jets one should 
choose fi ~ mz/3. A conservative approach is to vary the scale in an appropriate range and take the spread of the 
result as an estimate of the error due to higher order corrections. From fig. |2| we see that if we take /i in the range 
m.z/10 — mz the error due to the scale and a s in the determination of rhb{mz) would be of about 0.20 GeV . If scales 
as low as [i = 5 GeV are accepted the error increases to 0.23 GeV . Whether this error can be reduced by a clever 
choice of the scale or resummation of leading logs in y c or rb remains to be seen. 



The calculation presented in this paper has already been used by the DELPHI collaboration 16 to extract fhb(mz)- 
The preliminary result [we like to thank the DELPHI collaboration for allowing us to use these numbers here] they 
obtain is 

fh b (mz) = 2.85 ± 0.22 (stat) ± 0.20 (theo) ± 0.36 (fragmentation) GeV (9) 
which has to be compared with low energy determinations of the bottom quark mass. The last analysis of the T 



system using QCD sum rules |31 gives mb(fhb) = 4.13 ± 0.06 GeV which translates into fhb{mz) = 2.83 ± 0.10 GeV 
if one uses three loop renormalization group running and ct s (mz) = 0.118 ±0.003. On the other hand, the last lattice 
result gives |J2| mb{fhb) = 4.15 ± 0.20 GeV and fhb(mz) — 2.84 ± 0.21 GeV . Given the errors it is clear that central 
values agree so well just by chance. In particular the value in eq. (|^) has been extracted by using a different central 
value for a s . 

It is encouraging to see that this preliminary measurement is in full compatibility with low energy data, and, 
although for the moment it is not competitive with low energy measurements, it is good enough for testing the 
running of the bottom quark mass from /i = Mj, to /i = mz- the result for fhb(mz) in eq. (0) and the previous values 
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for mb{m\y) differ by more than 2.5 standard deviations. We believe that these results can be substantially improved 
with more experimental and theoretical work. 
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FIG. 2. Extracted value of mb(mz) if R^ixp — 0.96 as a function of the scale fi. We take a s (mz) = 0.118 (solid) and 
Aa s = 0.003 (dashed). 
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